home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_I0.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.0 KB  |  233 lines
  1. /* specfunc/bessel_I0.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_bessel.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "chebyshev.h"
  30. #include "cheb_eval.c"
  31.  
  32. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33.  
  34.  
  35. /* based on SLATEC besi0 */
  36.  
  37. /* chebyshev expansions
  38.  
  39.  series for bi0        on the interval  0.        to  9.00000d+00
  40.                     with weighted error   2.46e-18
  41.                      log weighted error  17.61
  42.                    significant figures required  17.90
  43.                     decimal places required  18.15
  44.  
  45.  series for ai0        on the interval  1.25000d-01 to  3.33333d-01
  46.                     with weighted error   7.87e-17
  47.                      log weighted error  16.10
  48.                    significant figures required  14.69
  49.                     decimal places required  16.76
  50.  
  51.  
  52.  series for ai02       on the interval  0.        to  1.25000d-01
  53.                     with weighted error   3.79e-17
  54.                      log weighted error  16.42
  55.                    significant figures required  14.86
  56.                     decimal places required  17.09
  57. */
  58.  
  59. static double bi0_data[12] = {
  60.   -.07660547252839144951,
  61.   1.92733795399380827000,
  62.    .22826445869203013390, 
  63.    .01304891466707290428,
  64.    .00043442709008164874,
  65.    .00000942265768600193,
  66.    .00000014340062895106,
  67.    .00000000161384906966,
  68.    .00000000001396650044,
  69.    .00000000000009579451,
  70.    .00000000000000053339,
  71.    .00000000000000000245
  72. };
  73. static cheb_series bi0_cs = {
  74.   bi0_data,
  75.   11,
  76.   -1, 1,
  77.   11
  78. };
  79.  
  80. static double ai0_data[21] = {
  81.    .07575994494023796, 
  82.    .00759138081082334,
  83.    .00041531313389237,
  84.    .00001070076463439,
  85.   -.00000790117997921,
  86.   -.00000078261435014,
  87.    .00000027838499429,
  88.    .00000000825247260,
  89.   -.00000001204463945,
  90.    .00000000155964859,
  91.    .00000000022925563,
  92.   -.00000000011916228,
  93.    .00000000001757854,
  94.    .00000000000112822,
  95.   -.00000000000114684,
  96.    .00000000000027155,
  97.   -.00000000000002415,
  98.   -.00000000000000608,
  99.    .00000000000000314,
  100.   -.00000000000000071,
  101.    .00000000000000007
  102. };
  103. static cheb_series ai0_cs = {
  104.   ai0_data,
  105.   20,
  106.   -1, 1,
  107.   13
  108. };
  109.  
  110. static double ai02_data[22] = {
  111.    .05449041101410882,
  112.    .00336911647825569,
  113.    .00006889758346918,
  114.    .00000289137052082,
  115.    .00000020489185893,
  116.    .00000002266668991,
  117.    .00000000339623203,
  118.    .00000000049406022,
  119.    .00000000001188914,
  120.   -.00000000003149915,
  121.   -.00000000001321580,
  122.   -.00000000000179419,
  123.    .00000000000071801,
  124.    .00000000000038529,
  125.    .00000000000001539,
  126.   -.00000000000004151,
  127.   -.00000000000000954,
  128.    .00000000000000382,
  129.    .00000000000000176,
  130.   -.00000000000000034,
  131.   -.00000000000000027,
  132.    .00000000000000003
  133. };
  134. static cheb_series ai02_cs = {
  135.   ai02_data,
  136.   21,
  137.   -1, 1,
  138.   11
  139. };
  140.  
  141.  
  142. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  143.  
  144. int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result)
  145. {
  146.   double y = fabs(x);
  147.  
  148.   /* CHECK_POINTER(result) */
  149.  
  150.   if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
  151.     result->val = 1.0 - y;
  152.     result->err = 0.5*y*y;
  153.     return GSL_SUCCESS;
  154.   }
  155.   else if(y <= 3.0) {
  156.     const double ey = exp(-y);
  157.     gsl_sf_result c;
  158.     cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
  159.     result->val = ey * (2.75 + c.val);
  160.     result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err;
  161.     return GSL_SUCCESS;
  162.   }
  163.   else if(y <= 8.0) {
  164.     const double sy = sqrt(y);
  165.     gsl_sf_result c;
  166.     cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c);
  167.     result->val  = (0.375 + c.val) / sy;
  168.     result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
  169.     result->err += c.err / sy;
  170.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  171.     return GSL_SUCCESS;
  172.   }
  173.   else {
  174.     const double sy = sqrt(y);
  175.     gsl_sf_result c;
  176.     cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c);
  177.     result->val = (0.375 + c.val) / sy;
  178.     result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
  179.     result->err += c.err / sy;
  180.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  181.     return GSL_SUCCESS;
  182.   }
  183. }
  184.  
  185.  
  186. int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result)
  187. {
  188.   double y = fabs(x);
  189.  
  190.   /* CHECK_POINTER(result) */
  191.  
  192.   if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
  193.     result->val = 1.0;
  194.     result->err = 0.5*y*y;
  195.     return GSL_SUCCESS;
  196.   }
  197.   else if(y <= 3.0) {
  198.     gsl_sf_result c;
  199.     cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
  200.     result->val  = 2.75 + c.val;
  201.     result->err  = GSL_DBL_EPSILON * (2.75 + fabs(c.val));
  202.     result->err += c.err;
  203.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  204.     return GSL_SUCCESS;
  205.   }
  206.   else if(y < GSL_LOG_DBL_MAX - 1.0) {
  207.     const double ey = exp(y);
  208.     gsl_sf_result b_scaled;
  209.     gsl_sf_bessel_I0_scaled_e(x, &b_scaled);
  210.     result->val  = ey * b_scaled.val;
  211.     result->err  = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val);
  212.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  213.     return GSL_SUCCESS;
  214.   }
  215.   else {
  216.     OVERFLOW_ERROR(result);
  217.   }
  218. }
  219.  
  220. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  221.  
  222. #include "eval.h"
  223.  
  224. double gsl_sf_bessel_I0_scaled(const double x)
  225. {
  226.   EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); ) 
  227. }
  228.  
  229. double gsl_sf_bessel_I0(const double x)
  230. {
  231.   EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); ) 
  232. }
  233.